/*
Copyright (c) [2019年5月1日] [吴超]
[MBT_studio] is licensed under Mulan PSL v2.
You can use this software according to the terms and conditions of the Mulan PSL v2.
You may obtain a copy of Mulan PSL v2 at:
		 http://license.coscl.org.cn/MulanPSL2
THIS SOFTWARE IS PROVIDED ON AN "AS IS" BASIS, WITHOUT WARRANTIES OF ANY KIND,
EITHER EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO NON-INFRINGEMENT,
MERCHANTABILITY OR FIT FOR A PARTICULAR PURPOSE.
See the Mulan PSL v2 for more details.
*/



#pragma once

//#define __CUDACC__
#include "graph_Adjacency.h"


#include <device_functions.h>
#include <device_launch_parameters.h>
#include <cuda_runtime.h>
#include <cuda.h>


#include <thrust/device_vector.h>
#include <thrust/sort.h>
#include <thrust/reduce.h>
#include <thrust/extrema.h>


//#include "几何图形/图形相交.h"
//#include <amp.h>
//#include <amp_math.h>

//#include <device_atomic_functions.hpp>
//#include <atomic>
#include <几何数据.h>
#include <Vector/向量运算.h>
#include <map>
#include <stdio.h>

#include <iostream>
using namespace std;




//texture<uint32, cudaTextureType2D, cudaReadModeNormalizedFloat> size_text;


/*template<typename T>
T**** f_create_3d_Array(uvec3 维度) {
	T* base;
	cudaError_t err = cudaMallocManaged(&base, 维度.x * 维度.y * 维度.z * sizeof(T));
	assert(err == cudaSuccess);
	T*** ary;
	err = cudaMallocManaged(&ary, (a + a * b + a * b * c) * sizeof(T*));
	assert(err == cudaSuccess);
	for (int i = 0; i < a; i++) {
		ary[i] = (T***)((ary + a) + i * b);
		for (int j = 0; j < b; j++) {
			ary[i][j] = (T**)((ary + a + a * b) + i * b * c + j * c);
			for (int k = 0; k < c; k++)
				ary[i][j][k] = base + ((i * b + j) * c + k) * d;
		}
	}
	return ary;
}*/


/*struct Lock {
	int* mutex;
	Lock(void) {
		cudaMalloc((void**)&mutex, sizeof(int));
		cudaMemset(mutex, 0, sizeof(int));
	}
	~Lock(void) {
		cudaFree(mutex);
	}
	__device__ void lock(void) {
		while (atomicCAS(mutex, 0, 1) != 0);
	}
	__device__ void unlock(void) {
		atomicExch(mutex, 0);
	}
};*/

#define DEF_线程数 1024

Inline S_计算维度 f_计算维度(uint32 v1num) {
	//uint32 最大粒子数 = DEF_Max(v1num, v2num);
	uint32 最大粒子数 = v1num;

	uint32 线程数 = DEF_线程数;
	//dim3 计算维度3D = { 1, 1, 1 };
	dim3 计算维度3D = { 1, 1, 1 };
	计算维度3D.y = 最大粒子数 / 线程数 + 1;
	if (计算维度3D.y > 线程数) {
		计算维度3D.y = 计算维度3D.y / 线程数 + 1;
		计算维度3D.x = 线程数;
	}
	uint32 计算维度2D = f_平方根(最大粒子数) + 1;


	uint32 内存用量3D = 计算维度3D.x * 计算维度3D.y * 线程数;
	uint32 内存用量2D = 计算维度2D * 计算维度2D;

	uint64 内存用量 = DEF_Max(内存用量3D, 内存用量2D);

	return { 计算维度3D , {计算维度2D, 计算维度2D}, 内存用量 };
}



__device__ inline uint32 f_CoordToHash(const float3& vert) {
	return uint32(vert.x + 511) | (uint32(vert.y + 511) << 10) | (uint32(vert.z + 511) << 20);
}
__device__ inline ushort3 f_CoordToVoxelCoord(const vec3& vert) {
	return { uint16(vert.x + 511) , uint16(vert.y + 511), uint16(vert.y + 511) };
}
__device__ inline ushort3 f_HashToCoord(const uint32& vert) {
	return { uint16((vert & 0x3ff) - 511) , uint16(((vert >> 10) & 0x3ff) + 511) | uint16(((vert << 20) & 0x3ff) + 511) };
}
__device__ inline uint32 f_CoordToIndex(const ushort3& coord, const ushort2& XZ维度) {
	return coord.z * XZ维度.x * XZ维度.y + coord.y * XZ维度.x + coord.x;
}



__global__ void f_KeyBuild(vec3* 顶点, uvec2* keys, float32 单元大小) {
	uint32 id = threadIdx.x + (blockDim.x * (gridDim.x * blockIdx.y + blockIdx.x));
	//uint32 id = blockIdx.x * blockDim.x * blockDim.y + threadIdx.y * blockDim.x + threadIdx.x;
	vec3 vert = 顶点[id];
	vert.x /= 单元大小;
	vert.y /= 单元大小;
	vert.z /= 单元大小;

	uint32 key = f_PackData3X32f(vert);
	//uint32 key = uint32(vert.x + 511) | (uint32(vert.y + 511) << 10) | (uint32(vert.z + 511) << 20);
	keys[threadIdx.x].x = threadIdx.x;
	keys[threadIdx.x].y = id;

}


__global__ void f_RadixSort(uint8 bit, const uint32 批量数, uvec2* 排序后对应键值, uvec2* 排序临时值A, uvec2* 排序临时值B, uvec2* count) {
	//uint32 id = threadIdx.x + (blockDim.x * (gridDim.x * blockIdx.y + blockIdx.x));
	uint32 id = (gridDim.x * blockIdx.y + blockIdx.x) * 批量数;
	uint32 id_x = (gridDim.x * blockIdx.y + blockIdx.x);


	uint32 桶1计数 = 0;
	uint32 桶2计数 = 0;

	uint32 bit_mask = (0x1 << bit);

	for (uint16 i = 0; i < 批量数; ++i) {
		uvec2 temp = 排序后对应键值[id + i];

		if ((temp.x & bit_mask)) {
			排序临时值A[id + 桶1计数] = temp;
			++桶1计数;

		}
		else {
			排序临时值B[id + 桶2计数] = temp;
			++桶2计数;
		}
	}

	count[id_x].x = 桶1计数;
	count[id_x].y = 桶2计数;

}

__global__ void f_RadixSort_Offset(const uint32 元素数量, uvec2* offset, uvec2* count) {
	uint32 桶1计数 = 0;
	uint32 桶2计数 = 0;

	桶1计数 = 0;
	桶2计数 = 0;

	for (uint32 i = 0; i < 元素数量; ++i) {
		offset[i].x = 桶1计数;
		桶1计数 += count[i].x;

		offset[i].y = 桶2计数;
		桶2计数 += count[i].y;
	}

}

__global__ void f_RadixSort_Merge(const uint32 批量数, const uint32 元素数量, uvec2* 排序后对应键值, uvec2* 排序临时值A, uvec2* 排序临时值B, uvec2* offset, uvec2* count) {
	uint32 桶1计数 = 0;
	uint32 桶2计数 = 0;

	uint32 id = (gridDim.x * blockIdx.y + blockIdx.x) * 批量数;
	uint32 id_x = (gridDim.x * blockIdx.y + blockIdx.x);


	桶1计数 = count[id_x].x;
	uint32 偏移 = offset[id_x].x;
	for (uint16 i = 0; i < 桶1计数; ++i) {
		//排序后对应键值[id + i] = 桶1[offset];
		排序后对应键值[i + 偏移] = 排序临时值A[id + i];
		//++offset;
	}

	桶2计数 = count[id_x].y;
	偏移 = offset[id_x].y + offset[元素数量 - 1].x + count[元素数量 - 1].x;
	for (uint16 i = 0; i < 桶2计数; ++i) {
		//排序后对应键值[id + 桶1计数 + i] = 桶2[offset];
		排序后对应键值[i + 偏移] = 排序临时值B[id + i];
		//++offset;
	}

}

















//邻接元素构建
__global__ void f_AdjacencyBuildVoxel(uvec2* 排序后对应键值, uvec2* offset, uint32* 体素索引, const ushort2 维度) {
	uint32 id = threadIdx.x + (blockDim.x * (gridDim.x * blockIdx.y + blockIdx.x));
	//uint32 id = (gridDim.x * blockIdx.y + blockIdx.x);
	if (排序后对应键值[id].x - 1 != 排序后对应键值[id].x) {
		auto coord = f_HashToCoord(排序后对应键值[id].x);
		体素索引[f_CoordToIndex(coord, 维度)] = id;
	}
}

__global__ void f_clearVoxel(uint32* 体素索引) {
	uint32 id = threadIdx.x + (blockDim.x * (gridDim.x * blockIdx.y + blockIdx.x));
	体素索引[id] = 0;
}



__device__ uint32 f_AdjacencyBuildIndex(const vec3 vert,
	const vec3* 邻接顶点,
	const uvec2* 排序后对应键值,
	const uint32* 体素索引,
	const uint32 voxelIndex[3][3],

	uint32* 紧凑键值,
	uint32 id,
	float32 粒子直径,
	uint8 邻接最大数量
) {
	uint32 当前邻接构建数量 = 0;

	for (uint8 x = 0; x < 3; ++x) {
		for (uint8 y = 0; y < 3; ++y) {
			auto 邻接顶点索引偏移 = 体素索引[voxelIndex[x][y]];

			uint32 开头坐标编码 = 排序后对应键值[邻接顶点索引偏移].x;

			for (uint32 i = 0; i < 邻接最大数量; ++i) {
				uvec2 邻接索引 = 排序后对应键值[邻接顶点索引偏移 + i];
				//uint32 邻接顶点索引 = 邻接索引.y;
				auto adjacencyVert = 邻接顶点[邻接索引.y];
				float32 len = vec_len(vec_sub(adjacencyVert, vert));

				if (len < 粒子直径) {
					//当前邻接紧凑键值[threadIdx.x][threadIdx.y] = 邻接顶点索引;
					紧凑键值[id * 邻接最大数量 + 当前邻接构建数量] = 邻接索引.y;
					++当前邻接构建数量;
				}


				if (当前邻接构建数量 >= 邻接最大数量) {
					return 当前邻接构建数量;
				}

				if (邻接索引.x != 开头坐标编码) {
					break;
				}
			}
		}
	}

	return 当前邻接构建数量;
}
__global__ void f_AdjacencyBuildElements(const uvec2* 排序后对应键值, const uint32* 体素索引,
	const vec3* 当前顶点, const vec3* 邻接顶点,
	uint32* 紧凑键值, uint8* 单元键值数量,
	const uint8 邻接最大数量, const float32 粒子直径, const ushort2 维度
) {


	//uint32 id = threadIdx.x + (blockDim.x * (gridDim.x * blockIdx.y + blockIdx.x));
	uint32 id = gridDim.x * blockIdx.y + blockIdx.x;
	//uint32 id = (gridDim.x * blockIdx.y + blockIdx.x);
	auto vert = 当前顶点[id];
	auto coord = f_CoordToVoxelCoord(vert);
	//auto coord	= f_HashToCoord(排序后对应键值[id].x);
	//auto 当前VertxIndex = 排序后对应键值[id].y;

	if (coord.x) {
		coord.x += threadIdx.x - 1;
	}
	if (coord.y) {
		coord.y += threadIdx.y - 1;
	}
	if (coord.z) {
		coord.z += threadIdx.z - 1;
	}


	__shared__ uint32 voxelIndex[3][3];
	//auto voxelIndex = f_CoordToIndex(coord, 维度);
	voxelIndex[threadIdx.x][threadIdx.y] = f_CoordToIndex(coord, 维度);
	__syncthreads();


	if (threadIdx.x == 0) {
		//uint32 当前坐标哈希 = 排序后对应键值[id].x;
		uint32 当前邻接构建数量 = f_AdjacencyBuildIndex(vert, 邻接顶点, 排序后对应键值, 体素索引, voxelIndex, 紧凑键值, id, 粒子直径, 邻接最大数量);
		单元键值数量[id] = DEF_Min(邻接最大数量, 当前邻接构建数量);

	}
	//单元键值数量[id] = 3;
}


template<typename T>
__device__ inline void f_Swap(T& a, T& b) {
	T temp = b;
	b = a;
	a = temp;
}

//冒泡排序
__global__ void f_BubbleSort(const uint32 数量, const uint32 维度, uint32* 排序键值, uint8* 是否结束) {
	uint32 num = 数量;

	uint32 idx_1 = threadIdx.x * 2 + 0 + (维度 * 维度 * blockIdx.y) + (维度 * blockIdx.x);
	uint32 idx_2 = threadIdx.x * 2 + 1 + (维度 * 维度 * blockIdx.y) + (维度 * blockIdx.x);
	uint32 idx_1B = threadIdx.x * 2 + 1 + (维度 * 维度 * blockIdx.y) + (维度 * blockIdx.x);
	uint32 idx_2B = threadIdx.x * 2 + 2 + (维度 * 维度 * blockIdx.y) + (维度 * blockIdx.x);

	for (uint32 i = 0; i < num; ++i) {
		排序键值[idx_1] = 0;
		if (排序键值[idx_1] > 排序键值[idx_2]) {
			f_Swap(排序键值[idx_1], 排序键值[idx_2]);
			是否结束[0] = 1;
		}
		//__syncthreads();

		if (排序键值[idx_1B] > 排序键值[idx_2B]) {
			f_Swap(排序键值[idx_1B], 排序键值[idx_2B]);
			是否结束[0] = 1;
		}
		//__syncthreads();
	}



}

/*__global__ void f_MergeSort(const uint32 数量, uint32* 排序键值A, uint32* 排序键值B) {
	uint32 id = (blockDim.x * blockDim.y + threadIdx.x) + (blockIdx.x * gridDim.y + blockIdx.y);
	uint32 tid = blockDim.x * blockDim.y + threadIdx.x;

	extern __shared__ uint32 shared_arr[];
	shared_arr[tid] = 排序键值A[id];
	__syncthreads();


	//__threadfence();
}*/



#define COUDA_CALL(x) {const cudaError_t a = (x);			\
	if (a != cudaSuccess) {									\
		cout<<"cuda err:"<<cudaGetErrorString(a)<<a<<endl;	\
		cudaDeviceReset();									\
		assert(0);											\
	}														\
}															\


static void f_构建体素索引(vec3* 邻接顶点, S_cu邻接体素索引& 体素索引, S_计算维度 维度, S_邻接点属性 属性) {
	//uint32 邻接计数构建维度 = f_平方根(num) + 1;

	//uvec2 *临时顶点键值	= 0;
	uvec2* 排序顶点键值A = 0;
	uvec2* 排序顶点键值B = 0;

	//COUDA_CALL(cudaMalloc((void**)&临时顶点键值, sizeof(uvec2) * 内存用量));
	COUDA_CALL(cudaMalloc((void**)&排序顶点键值A, sizeof(uvec2) * 维度.m_内存用量));
	COUDA_CALL(cudaMalloc((void**)&排序顶点键值B, sizeof(uvec2) * 维度.m_内存用量));

	uvec2* 排序组偏移 = 0;
	uvec2* 排序组数量 = 0;
	COUDA_CALL(cudaMalloc((void**)&排序组偏移, sizeof(uvec2) * 维度.m_计算维度3D.x * 维度.m_计算维度3D.y));
	COUDA_CALL(cudaMalloc((void**)&排序组数量, sizeof(uvec2) * 维度.m_计算维度3D.x * 维度.m_计算维度3D.y));


	f_KeyBuild << < 维度.m_计算维度3D, DEF_线程数 >> > (邻接顶点, 体素索引.m_排序索引, 属性.m_单元大小);
	for (uint8 bit = 0; bit < 32; ++bit) {
		f_RadixSort << < 维度.m_计算维度3D, 1 >> > (bit, 维度.m_计算维度3D.x, 体素索引.m_排序索引, 排序顶点键值A, 排序顶点键值B, 排序组数量);
		f_RadixSort_Offset << <1, 1 >> > (维度.m_计算维度3D.x * 维度.m_计算维度3D.y, 排序组偏移, 排序组数量);
		f_RadixSort_Merge << < 维度.m_计算维度3D, 1 >> > (维度.m_计算维度3D.x, 维度.m_计算维度3D.x * 维度.m_计算维度3D.y, 体素索引.m_排序索引, 排序顶点键值A, 排序顶点键值B, 排序组偏移, 排序组数量);
		cout << "GPU迭代:" << uint32(bit) << endl;
	}

	f_AdjacencyBuildVoxel << < 维度.m_计算维度3D, DEF_线程数 >> > (体素索引.m_排序索引, 排序组偏移, 体素索引.m_体素索引, { 512, 512 });

	//vector<uvec2> a(维度.m_内存用量);
	//cudaMemcpy(a.data(), 体素索引.m_排序索引, 维度.m_内存用量 * sizeof(uint32), cudaMemcpyKind::cudaMemcpyDeviceToHost);

	//cudaFree(临时顶点键值);
	cudaFree(排序顶点键值A);
	cudaFree(排序顶点键值B);
	cudaFree(排序组偏移);
	cudaFree(排序组数量);

}



void f_构建邻接索引(vec3* vert, vec3* 邻接顶点, S_cu邻接体素索引& voxel, S_cu邻接点索引组& 邻接索引, S_计算维度 dim, S_邻接点属性 属性) {

	f_AdjacencyBuildElements << < dim.m_计算维度2D, { 3,3 } >> > (voxel.m_排序索引, voxel.m_体素索引, vert, 邻接顶点, 邻接索引.m_邻接索引, 邻接索引.m_邻接数量, 属性.m_邻接最大数量, 属性.m_点半径 * 2, { 512, 512 });

}




void f_graph_分组顶点构建邻接信息(vec3* vert1, uint64 v1num, vec3* vert2, uint64 v2num, S_cu邻接点索引组& D邻接索引, S_cu邻接点索引组& S邻接索引, S_邻接点属性 属性) {
	/*uint32 最大粒子数 = DEF_Max(v1num, v2num);

	uint32 线程数 = 1024;
	dim3 计算维度3D = { 1, 1, 1 };
	计算维度3D.y = 最大粒子数 / 线程数 + 1;
	if (计算维度3D.y > 线程数) {
		计算维度3D.y = 计算维度3D.y / 线程数 + 1;
		计算维度3D.x = 线程数;
	}
	uint32 计算维度2D = f_平方根(最大粒子数) + 1;

	uint32 内存用量3D = 计算维度3D.x * 计算维度3D.y * 线程数;
	uint32 内存用量2D = 计算维度2D * 计算维度2D;

	uint64 内存用量 = DEF_Max(内存用量3D, 内存用量2D);*/

	S_计算维度 D动态粒子计算维度 = f_计算维度(v1num);
	S_计算维度 S动态粒子计算维度 = f_计算维度(v2num);



	vec3* 顶点1坐标 = 0;
	vec3* 顶点2坐标 = 0;

	uint32* 邻接索引1 = 0;
	uint8* 邻接数量1 = 0;
	uint32* 邻接索引2 = 0;
	uint8* 邻接数量2 = 0;

	uint32* 顶点体素 = 0;

	uvec2* 顶点1排序 = 0;
	uvec2* 顶点2排序 = 0;



	COUDA_CALL(cudaMalloc((void**)&顶点1坐标, sizeof(vec3) * D动态粒子计算维度.m_内存用量));
	COUDA_CALL(cudaMalloc((void**)&顶点2坐标, sizeof(vec3) * S动态粒子计算维度.m_内存用量));

	COUDA_CALL(cudaMalloc((void**)&邻接索引1, sizeof(uint32) * D动态粒子计算维度.m_内存用量 * 64));
	COUDA_CALL(cudaMalloc((void**)&邻接数量1, sizeof(uint8) * D动态粒子计算维度.m_内存用量));
	COUDA_CALL(cudaMalloc((void**)&邻接索引2, sizeof(uint32) * D动态粒子计算维度.m_内存用量 * 64));
	COUDA_CALL(cudaMalloc((void**)&邻接数量2, sizeof(uint8) * D动态粒子计算维度.m_内存用量));

	COUDA_CALL(cudaMalloc((void**)&顶点体素, sizeof(uint32) * 512 * 512 * 512));

	COUDA_CALL(cudaMalloc((void**)&顶点1排序, sizeof(uvec2) * D动态粒子计算维度.m_内存用量));
	COUDA_CALL(cudaMalloc((void**)&顶点2排序, sizeof(uvec2) * S动态粒子计算维度.m_内存用量));



	COUDA_CALL(cudaMemcpy(顶点1坐标, vert1, (sizeof(vec3) * v1num), cudaMemcpyKind::cudaMemcpyHostToDevice));
	COUDA_CALL(cudaMemcpy(顶点2坐标, vert2, (sizeof(vec3) * v2num), cudaMemcpyKind::cudaMemcpyHostToDevice));

	//COUDA_CALL(cudaMemcpy(邻接索引1, D邻接索引.m_邻接索引, (sizeof(uint32) * v1num * 属性.m_邻接最大数量), cudaMemcpyKind::cudaMemcpyHostToDevice));
	//COUDA_CALL(cudaMemcpy(邻接数量1, D邻接索引.m_邻接数量, (sizeof(uint8) * v1num), cudaMemcpyKind::cudaMemcpyHostToDevice));
	//COUDA_CALL(cudaMemcpy(邻接索引2, S邻接索引.m_邻接索引, (sizeof(uint32) * v1num * 属性.m_邻接最大数量), cudaMemcpyKind::cudaMemcpyHostToDevice));
	//COUDA_CALL(cudaMemcpy(邻接数量2, S邻接索引.m_邻接数量, (sizeof(uint8) * v1num), cudaMemcpyKind::cudaMemcpyHostToDevice));


	S_cu邻接体素索引	体素索引 = { 顶点1排序, 顶点体素 };
	S_cu邻接点索引组	邻接索引组 = { 邻接索引1, 邻接数量1 };
	{
		//动态粒子
		f_构建体素索引(顶点1坐标, 体素索引, D动态粒子计算维度, 属性);
		f_构建邻接索引(顶点1坐标, 顶点1坐标, 体素索引, 邻接索引组, D动态粒子计算维度, 属性);

		//vector<uint8> a(D动态粒子计算维度.m_内存用量);
		//cudaMemcpy((void*)a.data(), 邻接数量1, v1num * sizeof(uint8), cudaMemcpyKind::cudaMemcpyDeviceToHost);
	}

	f_clearVoxel << < {512, 512}, 512 >> > (顶点体素);

	{
		//静态粒子
		体素索引 = { 顶点2排序, 顶点体素 };
		邻接索引组 = { 邻接索引2, 邻接数量2 };

		f_构建体素索引(顶点2坐标, 体素索引, S动态粒子计算维度, 属性);
		f_构建邻接索引(顶点1坐标, 顶点2坐标, 体素索引, 邻接索引组, D动态粒子计算维度, 属性);


	}



	COUDA_CALL(cudaMemcpy(D邻接索引.m_邻接索引, 邻接索引1, v1num * sizeof(uint32) * 属性.m_邻接最大数量, cudaMemcpyKind::cudaMemcpyDeviceToHost));
	COUDA_CALL(cudaMemcpy(D邻接索引.m_邻接数量, 邻接数量1, v1num * sizeof(uint8), cudaMemcpyKind::cudaMemcpyDeviceToHost));
	COUDA_CALL(cudaMemcpy(S邻接索引.m_邻接索引, 邻接索引2, v1num * sizeof(uint32) * 属性.m_邻接最大数量, cudaMemcpyKind::cudaMemcpyDeviceToHost));
	COUDA_CALL(cudaMemcpy(S邻接索引.m_邻接数量, 邻接数量2, v1num * sizeof(uint8), cudaMemcpyKind::cudaMemcpyDeviceToHost));
	//f_graph_构建邻接信息(vert2, 0, v2num, 单元大小);


	cudaFree(顶点1坐标);
	cudaFree(顶点2坐标);

	cudaFree(邻接索引1);
	cudaFree(邻接数量1);
	cudaFree(邻接索引2);
	cudaFree(邻接数量2);

	cudaFree(顶点体素);

	cudaFree(顶点1排序);
	cudaFree(顶点2排序);
}








